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Abstract. The evolution of several physical and biological systems, ranging from 
neutron transport in multiplying media to epidemics or population dynamics, can be 
described in terms of branching exponential flights, a stochastic process which couples 
a Galton- Watson birth-death mechanism with random spatial displacements. Within 
this context, one is often called to assess the length £v that the process travels in a given 

■ region V of the phase space, or the number of visits ny to this same region. In this 

paper, we address this issue by resorting to the Feynman-Kac formalism, which allows 

^ ^ _ characterizing the full distribution of iy and ny and in particular deriving explicit 

■ moment formulas. Some other significant physical observables associated to iy and 
O , ny, such as the survival probability, are discussed as well, and results are illustrated 

by revisiting the classical example of the rod model in nuclear reactor physics. 
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1. Introduction 

Consider a single walker initially emitted from a point source at time tq = at position 
To, with velocity vq. Once emitted, the walker undergoes a sequence of displacements 
(at constant speed), separated by collisions with the surrounding medium. When the 
scattering centers encountered by the travelling particle are spatially uniform, the inter- 
collision lengths are exponentially distributed [H |2] , so that the displacements from r' 
to r in direction a; = v/|v| between any two collisions obey the probability density 



with V = |v| [3111]. The quantity a{r,v) represents the interaction rate per unit length 
and takes the name of total cross section: a{r, v) typically depends on the particle 
position and speed, and is proportional to the probability of particle-medium interaction 
along a straight line, carrying units of the inverse of a length [3j. At each collision, 
the incident particle disappears, and k particles (the descendants) are emitted with 
probability pfc(r, f), whose velocities are randomly redistributed in angle and intensity 
according to a given probability density Cfc(v' — )■ v|r), which in principle can vary as 
a function of the number k of descendants [3j. Each descendant will then behave as 
the mother particle, and undergo a new sequence of displacements and collisions, giving 
thus rise to a branched structure, as illustrated in Fig. [H As a particular case, when 
Pk = Sk,i we recover the well known Pearson random walk [H [2]. 

Branching random flights as described above lie at the heart of physical and 
biological modeling [5l |6], and are key to the description of neutron transport in 
multiplying media and nucleon cascades [7j, evolution of biological populations [8], 
diffusion of reproducing bacteria [9J, and mutation-propagation of genes [lO], just to 
name a few. For a detailed survey, ranging from the pioneering work by Galton 
and Watson on the extinction probability of birth-death processes to the recent 
developments, see, e.g., [5l[7]. 

A central question for random walks is to determine the occupation statistics of the 
stochastic paths in a given portion V of the phase space [HI [121 [131 [IH [13 [lEl [13 [TBI 
[T9l l20l [2T1 [22] . For exponential flights, the two natural observables of the system are 
the number ny of occurred visits to the volume V and the total length iy travelled in 
^ [31 [H [T71 [IHl [20]. In reactor physics, for instance, knowledge of iy allows assessing 
the neutron flux due to the chain reaction and hence the deposited power or the number 
of radiation-induced structural defects [3, [23]. In a model of epidemics outbreak, ny 
corresponds to the number of infections in a region as a function of the position of 
the initial infected person (as long as the number of infected people is small, so that 
nonlinear effects due to the depletion of the susceptibles can be neglected, and that 
spatial displacements can be described by a simple random walk [21]). The quantity 
ny occurs also in population genetics, where one might be interested in quantifying 
the number ny of mutations of a given kind V, starting from a single character, as a 
function of the number of generations (this is closely related to the Ewens' formula for 




(1) 



Figure 1. An illustration of branching exponential flights starting from a point source 
q and traversing a volume V in phase space. 



the mutation partition, when mutations are allowed to be recurrent [25]). The goal of 
this paper is to characterize the statistical properties of the number of visits ny and of 
the travelled lengths ^v foi' branching exponential flights. 

This paper is structured as follows. In Sec. |2]we first focus on the average quantities 
(£y) and {ny). Then, in Sec. [3] and H] we assess the full distribution of £y and ny, 
respectively, by resorting to the Feynman-Kac formalism. In particular, we show 
that recursive formulas for the higher moments of the travelled lengths and for the 
number of visits can be easily derived based on this approach. Some other significant 
physical observables associated to £y and ny are discussed in Sec. [H Then, in Sec. [6] 
we illustrate the proposed formalism on one-dimensional exponential flights, the so- 
called rod model, and support our findings with Monte Carlo simulations. Perspectives 



are finally presented in Sec. [71 Technical details are left to Appendix A , Appendix B 
and [Appendix C 



2. The average observables 

In the following, we introduce a few simplifying hypotheses, whose main advantage is 
to keep notation to a minimum, yet retaining the key physical mechanisms. Thus, 
we assume that displacements are performed at a constant speed v = Vq, i.e., that 
only the walker directions uj do change after collisions. We furthermore assume that 
the medium is spatially homogeneous, so that pk{r,v) and cr(r, t>) can be taken to 
be constant. Finally, we assume that the probability density of the directions for 
the outgoing particles is isotropic, independent of the number of emitted descendants. 
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Figure 2. Survival probabilities S^{x()) forpo = 0.2, pi = 0.3, andp2 = 0-5 (i^i = 1.3, 
and Lc — 3.9987). Blue crosses and red circles: S^{xo) and S^(xo), respectively, 
with xq = 1.75 and i = 2 (A < 0). Green squares and black triangles: S^{xo) 
and S^{xq), respectively, with xq = 3.75 and L = 6 (A > 0). Solid lines are numerical 
integrals of Eq.[34l symbols Monte Carlo simulations with 10^ histories. Dashed curves: 
asymptotic survival probabilities Soo from Eq. [36l 



namely, 

CkiV -> v|r) = C{uj' ^ a;|r) = (2) 

where the normahzation factor = 2n'^^^ /T{d/ 2) is the surface of the unit sphere in 
dimension d. 

Branching exponential flights, as defined above, are a Markovian stochastic process 
that can be observed both as a function of time r and discrete generations n (this 
latter case corresponds to recording the particle position and direction at collision 
events only). Markovianity is granted by the fact that displacements between collisions 
are exponentially distributed [31 SI [7], and implies that knowledge of the phase space 
variables r, uj at time r or generation n is sufficient to determine the future evolution 
of the walker 0. 

To begin with, we address first the average physical observables of branching 
exponential fiights, which in most cases provide a reasonable first-order estimate of 
the system evolution. 

2.1. The average total travelled length {iy) 

Let A^(r, uj, r|ro, u^o) be the average number of particles that at time r are found in the 
phase space element drduj around r, lj, starting with a single particle emitted at tq in 
direction cjq, i.e., A^(r, a;, 0|ro, <-i;o) = 1 = — ro)^(^ — '^o)- Consider a displacement 

X Conversely, knowledge of the position r alone does not ensure Markovianity, as would instead be the 
case for branching Brownian motion [26[ [27l [28l [29] . 
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Figure 3. Survival probabilities S^{xo) forpo = 0.2, pi = 0.3, andp2 = 0.5 (i/i = 1.3, 
and Lc = 3.9987). Blue crosses and red circles: S'+(a;o) and S'^(xo), respectively, 
with a;o — 1-75 and i = 2 (A < 0). Green squares and black triangles: S^{xo) 
and S~{xo), respectively, with Xq = 3.75 and L = 6 (A > 0). Solid lines are numerical 
integrals of Eq.[35l symbols Monte Carlo simulations with 10^ histories. Dashed curves: 
asymptotic survival probabilities from Eq. [36l 



in a small time dr along a line oriented as iv: the time variation of iV(r, a;, r|ro, o^o) 
reads 

d f 

—N = -vaN + vi j du'vaC{ui' ^ u\y)N{t,uj' ,t\to,ujo), (3) 

the quantity ui = J2k ^Pk being the average number of secondary particles emitted per 
collision event. This equation can be understood as a mass balance: in dr, N decreases 
because of particles that at a time rate va interact, and thus change direction, and 
increases because of particles that, travelling in another direction a;', have a collision, are 
multiplied by a factor i/i, and change their direction to a; [23]. By using = ^+vu:-Vr, 



dr dr 

we finally have 

which is a Boltzmann-like conservation equation for the average particle density in 
phase space. Boundary conditions on A^(r, lj, r|ro, o^o) depend on the specific problem 
under analysis. Actually, instead of it is often common to introduce the quantity 
(f) = Nv, which takes the name of particle flux [23]. The stationary behavior of 
the particle density is provided by integrating over time, and for the stationary flux 
0(r, Li;|ro, ujq) = dT(f){r, uj, r|ro, ujq) we get in particular 

<-j ■ Vr0 + cr0 = i/i J duj'aC{uJ>' uj\v)(f){v,Lj'\ro,u;o) + q. (5) 

Eq. [5] can be recast in the more compact formula 

£0(r,a;|ro,a;o) = -g, (6) 
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Figure 4. Asymptotic survival probabilities S^{xo) for po = 0.2, pi = 0.3, and 
P2 = 0.5 (z/i = 1.3, and Lc = 3.9987). Blue circles and red squares: S'^{xo) and 
S^{xo), respectively, with L = 6 (X > 0). Solid lines are numerical integrals from 
Eq. [36l symbols Monte Carlo simulations with 10^ histories. 



where C = —uj ■ Vr — ct + z^i J duj'aC{uj' — ?■ a;|r) takes the name of (forward) transport 
operator [3l |23]. The quantity 0(r, c^lro, uJo) can be intepreted as the stationary density 
of the total length travelled by the particles in the phase space element drdu: around 
r, u}-. hence, the average travelled length in a given volume V of phase space will be 
given by 

)(ro,a;o) = y dr j c?a;\/(r, a;)</)(r, a;|ro, a;o), (7) 

where ^(r, (jj) denotes the marker function of the phase space volume V , i.e., V{y, uj) = 1 
when r, u; belong to V, and V{v,(jj) = elsewhere. 

2.2. The average total number of visits (riy) 

If generations are considered instead of time, a different mass balance equation for 
exponential flights can be established. Let '?/'„(r, ci;|ro, c*^o) be the average number of 
particles that enter a collision at r, having direction uj, at the n-th generation. Then, 
the following recursive formula can be established 



J dr' j duj'T{r' ^ r\uj)C{uj' ^ uj\r')^n{r',uj'\ro,uJo), 



with the initial condition ^/'i(r, a;|ro, tt'o) = T(ro — t- r|ct;o) [23]. The term ipi 
represents the average particle number entering a collision at the first generation (the 
so-called uncollided density). The stationary behavior of '?/'n(r, <^|ro, <^o) is obtained by 
summing over all generations: as customary, we define the collision density as being 
ip{r,uj\rQ,ujQ) = Y^'^^^ipn{T',<^\^o,<^o) [23], and we thus get the integral equation 



= ui 



j dr' j duj'T{r' r\uj)C{uj' Lc;|r )^(r', Lj'|ro, ^^o) + V^i- (9) 



Figure 5. Average and second moment of ly^ for xq = 1.75 and L = 2, with po — 0-2, 
Pi = 0.3, "Pi = 0.5 (i^i = 1.3 and L < Lc). Blue crosses: {iy)^ {xo); red circles: 
(^y)^(xo). Green squares: (€y)^(a;o); black triangles: {ty)^{xo). Solid lines are 
numerical integrals from Eq. 1161 symbols Monte Carlo simulations with 10^ histories. 
Dashed lines are the asymptotic limits in Eq. [20l 



The quantity '?/'(r, a;|ro, Wq) physically represents the stationary density of the number 
of particles entering a collision at r, uj: then, the average number of visits to a given 
volume V of phase space will be given by 



From (nv)(ro, ujq) and (£y)(ro, ujq) being two average observables of the same stochastic 
process, and thus closely related to each other, one can easily imagine that the 
quantities ip and cj) must be intimately connected as well. Actually, it can be shown 
that ?/'(r, ci;|ro, cl»o) = o"0(r, a;|ro, cjq) [21 [23]: this provides the relation between the 
stationary densities A^, ip and 0, and implies in particular that Eq. [9]can be equivalently 
recast into Eq. [5l by setting ip = acj). As a consequence, we have also (ny)(ro, cjq) = 
f dr J dujV{r,u})a(l){r,uj\ro,u)Q). 

The approach proposed in this Section so to assess the behavior of the average 
quantities (iy) and (ny) can not be straightforwardly extended to higher moments 
(which are often necessary to quantify the statistical fluctuations around the average), 
nor to other observables. In the following, we show that this difficulty can be overcome 
by resorting to the Feynman-Kac formalism, which allows characterizing the full 
distribution of the variables iy and ny- 




(10) 



8 




Figure 6. Average and second moment of ny, for xo = 1.75 and L ^ 2, with po — 0.2, 
Pi = 0.3, p2 — 0.5 {i^i — 1.3 and L < Lc). Blue crosses: (fiy)^(a;o); red circles: 
("v)n(2^o)- Green squares: (ny) + (xo); black triangles: (ny)~(xo). Solid lines are 
numerical integrals from Eq. 1261 symbols Monte Carlo simulations with 10^ histories. 
Dashed lines are the asymptotic limits in Eq. [29l 



3. Total travelled length in V 

We formally define the total length travelled by a branching exponential flight in 
a given volume V of the phase space, when observed up to a time t, as 



it) 



[ V{r',u;')vdt', (11) 
Jo 

where the integral is intended over all the branching paths of a single realization up 
to time t. The quantity iyit) is clearly a stochastic variable, which depends on the 
realizations of the underlying process, as well as on the initial conditions. Instead 
of studying the probability density function P((£y |ro, <^o), it is more convenient to 
introduce the associated moment generating function 

gi(s|ro,a;o) = (e-^'-(*))(ro,a;o), (12) 

where s is the transformed variable with respect to iy- Basically, Qt{s\rQ,ujQ) can 
be interpreted as the Laplace transform of Pt{iv\^o,<^o)- We derive then a backward 
equation for (s|ro, uJq) by closely following the approach originally proposed by Kac for 
Brownian motion [30] , based on Feynman path integrals As detailed in Appendix A 



the resulting backward Feynman-Kac equation relates the generating function Qt of the 
travelled length iy to the generating function G[z] = 'Yl,kPkZ^ offspring number 

§ The Feynman-Kac formalism more generally applies to continuous-time Markov processes (see, 
e.g., [3TJ |32l |33l [34]), and has recently been extended to non-Markovian walks [35 l [36 l [37 l l38]. 
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k, and reads 

^^Qt = ojo ■ VroQt - crQt - sV{ro,ujo)Qt + cxG[C*{Qt}], (13) 
where C*{Qt} is a shorthand for the direct ion- averaged Qt 

C*{Qt} = J C*{iv', ^ a;o|ro)Qt(s|ro,a;'o)da;'o, (14) 

C* being the adjoint probabihty density with respect to C. Equation [13] is completed 
by the initial condition (5o(s|ro,Wo) = 1 and by the appropriate boundary conditions, 
which depend on the problem at hand. 

3.1. Moment equations 

Equation [13] is a partial differential equation with a nonlinear integral term, for which 
explicit solutions are hardly available. Moreover, one would still need to invert the 
solution Qt so to obtain the probability density of iy in the direct space. A somewhat 
simpler approach consists in deriving the corresponding moment equation^: by the 
definition of Qt, the moments of the travelled length can be obtained from 

(£™),(ro,a;o) = (-l)"^Q*(s|ro, a;o)|.=o. (15) 

By taking the m-th derivative of Eq. [13] and resorting to the Faa di Bruno's formula for 
multiple derivatives of composite functions [3H], we get the following recursive formula 
for the moments of the trace length 

1 fi 

^ i=2 
for m > 1, where 



(17) 



is the (backward) transport operator adjoint to C [23]. Here Bm,j [zi\ = 
Bm,j [zi, Z2, - ■ ■ , Zm-j+i] are the Bell's polynomials [39], and uj = {k{k—l)...{k—j+l)) are 
the falling factorial moments of the descendant number, with z/q = 1. Bell polynomials 
commonly appear in connection with the combinatorics of branched structures [39]: this 
might give a hint about their role in Eq. [161 which relates the moments {iy)t of the 
travelled length to the moments uj of the descendant number. The recurrence is initiated 
with the conditions {iv)t = 1 (from normalization), and {iv)o = 0- Observe that {iv)t 
depends only on z/i, {iy)t on z/i and z/2, and so on. 

II For Brownian motion without branching, a similar approach was proposed by Kac |33| and later 
extended, e.g., in |141 1151 116] . Similarly, the moments of exponential flights without branching are 
discussed, e.g., in [HI [H [H [501 [H] . 

^ The first few polynomials read: Sq.o = 1; 'Si^ilzi] — zi; B2,i[zi, Z2] — Z2,B2,2[zi, Z2] = 
z1;B3^i[zi,Z2,Z3] = Z3,;B3,2[zi, 22,23] = 3ziZ2,'B3,3[zi, 22,23] = zf;.... 
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3.2. Stationary behavior 

Most often, the observation time t is much longer than the characteristic time scale of 
the system dynamics, which means that trajectories are followed up to t — oo. In this 
case, the time derivative in Eq. [16] vanishes, provided that the moment {('■y)t does not 
diverge when t — )■ oo. We therefore get a recursive formula for the stationary moments 
=lim,^oo(0*. namely, 

r(£^)(ro,a;o) = -[/,„_i(ro, a;o), (18) 

where 

m 

f/™_i(ro,a;o) =mV(ro,a;o)(£r') + ^5Zi^j-^™,.- [C*{{^v)}] (19) 

i=2 

can be interpreted as a (known) source term that depends at most on the moments of 
order m — 1. Now, from C* being the adjoint operator with respect to C, if one can solve 
Eq. [6] for a point source q = 6{r — ro)6{uj — ujq) and obtain the corresponding stationary 
flux 0, then Eq. [18] can be explicitly inverted, and gives 

(f;?)(ro,a;o) = j dr j rfa;0(r, u;|ro, a;o)t/„_i(r, a;), (20) 

which means that the stationary moments of the travelled length can be obtained by 
convoluting the stationary flux with the source term Um-i- As a particular case, when 
Pk>2 = we reobtain the simpler recursive formula derived in [191 122] for non-branching 
exponential flights. 

Finally, for the average length travelled in V , i.e., m = 1, we recover the formula 
{iy) = J dr J cia;l^(r,a;)0(r, a;|ro,a;o), since Uo{r,uj) = V{r,u!). 



4. Total number of visits to V 

We address then the statistical properties of the total number of visits nv{n) performed 
by a branching exponential flight in a given volume V of the phase space, when observed 
up to the n-th generatiocEl- We formally define 

^yH = I]nr„a;,), (21) 

i 

where the sum is intended over all the points visited by the branching path up to entering 
the n-th generation. We adopt the convention that the source is not taken into account. 
The quantity nv{n) is again a stochastic variable depending on the realizations of the 
underlying process and on the initial conditions. Similarly as done for iy, instead of 
studying the probability P„(ny |ro, uJq), it is more convenient to introduce the associated 
moment generating function 

Q„(n|ro,a;o) = (e-™-("))(ro,a;o), (22) 

+ The analysis of discrete (isotropic) branching walks can be extended to arbitrary flight length densities 
T, as shown in [401 141j. by resorting to an integral formulation. Here we stuck however to exponential 
flights, which imposes T as given in Eq. [TJ 
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where u is the transformed variable with respect to ny. As shown in Appendix B 



the backward discrete Feynman-Kac equation for Q„(M|ro,ajo) relates the generating 
function Qn of the number of visits ny to the generating function G of the offspring 
number k, and reads 

-^o-Vrogn+i(M|ro,a;o)+agn+i(w|ro,a;o) = ae-"^(^»''^°)G [C*{Qn}] , (23) 

with the initial condition 

gi(M|ro,a;o) = j e-^'^^'^^'^'^T* {r, ^ ro^dr,, (24) 

and the appropriate boundary conditions. We have used the shorthand C*{Qn} = 
J C*{uJq — )■ ujQ\Ti)Qn{u\ri,ujQ)duj'Q. As a particular case, when particles can not move 
(i.e., a — )■ oo), spatial dependences can be neglected and from Eq. [23] one recovers 
the counting statistics of a simple Galton- Watson process observed up to the n-th 
generation [5]. 



Moment equations 

Equation [23] is a nonlinear integro-differential and finite differences equation. Similarly 
as in the case of iy, the analysis of the distribution of ny can be simplified by deriving 
the corresponding moment equations: by the definition of Qn, the moments of the 
number of visits can be obtained from 

{n^Uro^LJo) = i-ir-^Qniu\ro,u;o)\u=o- (25) 

Then, by taking the m-th derivative of Eq. [231 we get the following recursive formula 
for the moments of number of visits 

m 

- LJQ ■ Vro(^™)„+l + (r{n'^)n+l = ^Y1 ^^^"^'^ [C*{{K)n}] + 

3=1 

m y \ m—k 

(r)^(^0''^0) E ^3^rn-k, [C*{{n^y)n]] , (26) 
k=l ^ ^ i=0 

for m > 1. Equation [26] relates the moments {ny)n of the number of visits to the 
moments Uj of the descendant number. Observe that {ny)n depends only on ui, {ny)n 
on Vi and V2i and so on. The recurrence is initiated with the conditions {ny)n = 1 (from 
normalization), and {ny)i = J V{ri,ujQ)T*{ri — )■ ro|u;o)(iri. 



4-2. Stationary behavior 

Most often, one considers trajectories that are followed up to n — )• oo, provided that the 
moment {ny)n does not diverge. We therefore get a recursive formula for the stationary 
moments (uy) = lim„_j.oo(^™)n, namely, 

C*{n^) {ro,LJo) = -ff^_i(ro,a;o), (27) 
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where 

m /■ \ m—k m 

H^_, = (7)^(^0' '^o) 5Z [C*{«)}]+a5^z/,S„,, [C*{{n^y)}] .(28) 

k=l ^ ^ i=0 i=2 

is a source term, and we have singled out the term of order m in the Bell 
polynomials. The quantity Hm-i is closely related to Um-i, and the contribution 
'^Y1]L2 ^j^rnj [C*{{ny)}] is common to both. When m = 1, we have Hq = alio. 

As done for the moments of travelled lengths, Eq. [27] can be explicitly inverted in 
terms of the corresponding stationary flux (p, and gives 

(n^)(ro,a;o) = y dv j rfa;0(r, a;|ro, a;o)i^m-i(r, a;), (29) 

which means that the stationary moments of the number of visits can be obtained by 
convoluting the stationary flux with the source term Hm-i- 

In particular, for the average number of visits to (m = 1) we recover the formula 
{riy) = J dr J dujV{rQ,u)o)a(f){r,uj\ro,u!o), since Ho{r,uj) = (rV{ro,uJo). The close 
relation between (riy) and (iy) carries over to moments of any order m > 1 (via Eq. [27|) . 
although the simple proportionality that holds true for the average does not apply to 
higher moments. 



5. Other physical observables 

The Feynman-Kac approach proposed in the previous Sections allows fully characterizing 
the distribution of the travelled length and of the number of visits. In the following, 
we show that a number of other interesting features of the process can be assessed by 
relying upon the same formalism, with minimal modifications, and we discuss some 
significant examples. 



5.1. Probability of never visiting a region V 

For instance, one might be interested in determining the probability -R„(ro,u^o) that 
a branching exponential flight coming from a point source at rQ,uJo never collides in 
a given domain V, up to generation n. This is intimately related to the well-known 
gambler's ruin problem [HE]. It follows that Rn{ro, loq) = Pniny = 0|ro, uJq)- Now, if we 
introduce the probability generating function (■u"^*^"))(ro, cjq) = J2i u^Pn{nv = i\ro, (^o), 
then by construction _R„(ro,a;o) = (M"'^''"^)(ro, a^o)U=o- By comparing (M"^*^"))(ro, cJq) 
to (5n(w|ro5 <-^o) = (e~""^''"^)(ro, cjo), it is apparent that /2„(i'05<^o) satisfies 

- ojo ■ Vro^n+i(ro,'^o) + cf Rn+i{r Q , o) = aV{ro,u}o)G [C* , (30) 

where V{ro,u!Q) = 1 — V{rQ,ujo). As for the initial conditions, we have Ri{rQ,ujQ) = 
J V{ri,ujQ)T*{ri — ro\uJo)dri. When n — )■ oo, we get the stationary probability 
equation 

- LJo ■ Vro/2(ro, + ai?(ro, u^o) = crVivo, u;o)G [C* {R}] , (31) 
where we have set -R(ro, ujq) = lim„_j.oo R 
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5.2. Survival probability 

The quantities iy and ny, as defined above, are cumulative, i.e., their knowledge 
integrates the whole history of the walk, from the source to the measure time, 
or generation. Sometimes, it is necessary to provide information about local 
(instantaneous) countings, namely the number mv{t) and my(n) of particles in 
the volume V when an observation is performed at time t, or at generation n, 



respectively (see, e.g., jTj). As shown in Appendix A the probability generating function 
Wt{s\ro,uJo) = (s"^^W) satisfies 

-^Wt{s\ro,oJo) = uJo-VroWMro,oJo)-aWt{s\ro,u:o)+cxG[C*{Wt}], (32) 
V at 

with initial condition M/o('S|ro, c<^o) = 3^^^°'^^°^ Equation [32] is known in reactor physics 



as the Pal-Bell equation |l2l |l3l [7]. Furthermore, as detailed in Appendix B the 
probability generating function H/'„(M|ro, cjq) = {u^^^"^^) satisfies 

- LJo ■ VroW„+i(M|ro, a;o) + aWn+i{u\vo, u^o) = aG , (33) 

with initial condition Wi{u\YQ,ijjQ) = J u^^^'^''^°^T*{ri — )■ ro|cL>o)(iri. 

Assume now that V is bounded, i.e., particles are lost upon leaving the volume: 
then, we might want to assess the survival probability at time t or generation n, due to 
the interplay between the branching mechanism and the spatial leakages. As particles 
can not re-enter V after crossing the boundaries, if my{t) = 0, then also my(t') = for 
t' > t (i.e., the process goes to extinction), and the same holds for mv{n). Hence, by 
definition, the probability of having zero particles in the volume at a time t is given 
by Wt{s = 0|ro,<-Jo), which equivalently yields the probability that extinction is reached 
for times smaller than t, since V is bounded [7]. We define then the survival probability 
as S'f(ro, ujq) = 1 — Wt{s = 0|ro, ujq), which by direct substitution in Eq. [32] satisfies 
1 d 

-^5'j(ro, ujq) = ujQ ■ Vro5't(ro, ujq) - (rSt{ro, ujq) - aF[C*{St}], (34) 

where we set F[z] = J2T=i^kz'' , with = {—l)''i'k/k\ [43j. At the boundaries, St must 
vanish when ojq is directed towards the exterior of V. The probability of having zero 
particles in the volume V at generation n is similarly given by Wn{u = 0|ro,ci;o), which 
therefore yields the probability that extinction is reached for generations smaller than n. 
We define then the associated survival probability as ^^(ro, ujq) = 1 — Wn{s = 0|ro, uJq), 
which by direct substitution in Eq. [33] satisfies 

<^o ■ Vro5'„+i(ro,a;o) - aSn+i{ro,ujo) = aF[C*{S„}] , (35) 

where S'„ must again vanish at the boundaries when ujq is directed towards the exterior of 
V. Finally, by either taking the limit S = limt_j.oo 5'((ro, ujq) or S = lim„_s.oo >S'„(ro, loq), 
respectively, the probability of ultimate survival ^(ro, uJq) satisfies 

ujo -"^roS -crS = aF[C*{S}]. (36) 

Observe that in principle S' = is always a solution to Eq. [361 which would imply 
almost sure extinction, i.e., a vanishing probability that infinitely long branching chains 
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exist in V. From the Galton- Watson theory, we know that when z/i < 1 the branching 
process goes to extinction even in the absence of spatial leakages, hence S = 0. However, 
when z/i > 1 the branching process would grow indefinitely, and it may happen that 
the particle loss due to finite geometry is not sufficient to compensate the population 
growth. In this case, the solution S = would become unstable, and St (or Sn) would 
converge towards a nontrivial survival probability 5* = 6*00 > 0. The stability analysis of 
the solution S = can be carried out by introducing a small perturbation, for instance 
in the form S ~ eX{t)Y{ro,ujo), the amplitude e > being a small positive constant, 
with X{t) > and Y{rQ,u:Q) > 0. Now, if we inject S into Eq. [Ml and take the limit 
e — 0, we obtain an equation for the perturbation amplitude 

1 1 dXjt) ^ ujQ ■ V.^y (rp, ojq) - (rp, ujq) + aiy,C*{Y} 

vX{t) dt Y . l-^O 

where at the numerator of the right hand side we recognize the adjoint operator £*. 
From the separation of the variables, Eq. [371 shows that the evolution of the perturbation 
amplitude with respect to time is determined by the ratio A = C*Y/Y , hence by the 
eigenvalue equation C*Y = XY. The spectrum of the eigenvalues of £* depends on the 
geometry of V and on the boundary conditions. If all eigenvalues A are negative, the 
amplitude of the small perturbation will shrink in time, so that eventually S = 0; if on 
the contrary at least one eigenvalue is positive, then the small perturbation will grow 
in time, which means that 5* = is unstable, and eventually 5* = Soo- For a given 
branching process, the crossover between these two regimes depends on the size and 
shape of the volume V, and generally speaking one would expect that S = Soo for a 
volume size larger than some critical value Vc [23l [7j, which is attained when A = 0. 
When ui < 1, Vc = 00. The stability analysis of S* ~ eX„F(rp, cjp) in Eq. [35] leads to 

^ = pc-w ^ 

Xn a;p ■ Vro>^(rp,a;p) - crF(rp,a;p) 

The evolution of the perturbation amplitude is determined then by the eigenvalue 
equation cjp ■ Vro^(rO)<^o) ~ cr^(ro!<^o) + ^C*{Y} = 0: when /3 < 1 the pertubation 
shrinks and when /3 > 1 the pertubation grows, the crossover occurring for a critical 
volume Vc such that (3 = 1. 

6. The rod model 

In order to illustrate the previous results, we revisit here a relevant example inspired 
by reactor physics. Neutrons in a multiplying medium undergo branching exponential 
ffights, where conceptually radiative capture represents absorption (po), scattering 
corresponds to pi and fission to pk>2'- in realistic situations, the cross sections and all 
the other physical parameters depend on energy and position, and angular distributions 
are often mildy or even strongly anisotropic [23]. To simplify the matter, we assume 
that cross sections are constant, scattering is isotropic and particles travel at an average 
constant speed v = 1. Moreover, we address a one- dimensional configuration, namely 
the interval [0, L], where only two angular directions (forward or backward) are allowed: 
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all physical quantities will be then denoted with a superscript ± according to whether 
(jj is taken along the x-axis (+), or in the opposite direction (— ). Despite these 
many simplifications, the so-called rod model yet captures the key features of neutron 
transport, and has been widely adopted in neutronics O HI]. 

To begin with, we impose leakage boundary conditions at x = and x = L 
and study the survival probabilities. By setting Y~^{xo) = Y{xo,Uo = +) and 
Y^{xo) = Y{xo,uJo = — ), the eigenvalue equation C*Y = XY gives rise to a system 
of two coupled linear differential equations 

^ -F+(xo) - ar+(xo) + ^ (r+(xo) + F-(xo)) = Ar+(xo), 



dxo 2 

-Y-{xo) - aY-{xo) + ^ (r+(xo) + ^-(xo)) = XY'^xo), (39) 



dxo 2 
where we have used C*{Y} = (F+(xo) + F^(xo))/2, thanks to isotropy. The general 
integrals Y~^{xo) = Yi(xo;ki,K2) and Y~{xo) = Y2{xo; ki, K2) of Eq. [39] can be easily 
obtained as a combination of exponential functions, up to two integration constants, 
say Ki and K2, that are to be imposed by boundary conditions. Due to the leakages, 
at the boundaries we have Y^{L) = and Y^{0) = 0. The solutions Y^ = Y^ = 
clearly satisfy the system in Eq. |39l together with boundary conditions. Searching for 
nontrivial solutions leads then to solving 

1^ F2(0;m,«:2) =0 J ^ ' 

with respect to the basis of ni and k,2- This in turn gives an implicit equation for the 
eigenvalues A as a function of the system parameters L, cr, and z/i, namely. 



/ X (Act - y) ^^"^^ f^o-A/MAT^^^x 

cosh (Lcta/ Act (Act - 1^1) ) H , = 

V / a/A_(A_ — 7/1) 



0, (41) 



a/Act (Act - l^l) 

where we have set the dimensionless variable Act = A/a + 1 . 

When A = A(L,cr, z/i) < 0, then S = and the neutron chain reaction will die 
out because of leakages and possibly absorptions; when on the contrary A > 0, then 
S = Soo{xo), and the chain reaction will diverge, as particles born from fission are not 
sufficiently compensated by leakages and absorptions; the crossover between these two 
regimes is reached for A = 0, which therefore defines the portion of the parameter space 
(L, 0", Ui) for which the branching process will attain ultimate extinction. When ui < 1, 
all A stay negative. When ui > 1, imposing A = in Eq. |4T] yields the explicit relation 

,-1 



L. = = ^A^^ , (42) 

cr - 1 

where Lc is the maximum system size such that for L > Lc there exists a finite survival 
probability that the neutron population will grow indefinitely. By resorting to the same 
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arguments, the analysis of the eigenvalue equation ujq ■ Vro^(ro)^o) ^ <^^(i"07'^o) + 



^C*{Y} = leads to 



(/3-f)sinh (La^l - ^) 



cosh ( - + '^,(;-;: ^ - (^3) 



which closely resembles Eq. SH Imposing /3 = 1 not surprisingly yields again Eq. US 
In a previous paper we derived Eq. |12] based on the analysis of the moment {uy) as a 
function of the system parameters [H]. 

The survival probabilities in Eq. [M] and [35] can be integrated numerically: in 
Figs. E] and [3] we compare the resulting curves (as a function of time or generations, 
respectively) with Monte Carlo simulations for different configurations. In particular, 
once the probabilities pk have been chosen, by varying the rod size L it is possible to 
impose L < Lc or L > L^. In the former case, the survival probabilities converge to zero 
independent of the starting direction, as expected, whereas in the latter the survival 
probabilities saturate to an asymptotic value 5*00 that depends on the starting point as 
well as on the initial direction of the walker. Observe that S'^ up to a time of the order 
of r+ ^ |L — xq\/v does not feel the effects of the boundaries, yet, and the same holds 
for Si up to a time of the order of r~ ~ xq/v. Therefore, we expect 5*^^ ^ up to 
min(r+,r^). For the configuration where L > Lc, the asymptotic survival probability 
as a function of starting point and initial direction is displayed in Fig. HI 

We examine then the statistics of the travelled length iy and number of visits 
ny via the moment equations [16] and [26] respectively. In particular, we consider here 
the average and second moment, which in most cases are sufficient to characterize the 
typical behaviour of the stochastic variables and their dispersion. The equations are 



given in Appendix C for a single neutron initially emitted at xq. We consider the same 



geometrical configuration as above, and take V{x,ijj) = V{x) for x G [0,L], i.e., we 
measure lengths and collisions independent of the local direction of the particles. As 
for boundary conditions, leakages impose (£y)^(L) = and (^y)r(O) = 0; similarly, 
(ny)^(L) = and {ny)~{0) = 0. For L > Lc the moments would diverge as time and 
generations increase. Therefore, we choose a configuration such that L < Lc and in 
Fig. [5] and [6] we plot the evolution of the moments of the travelled length and number 
of visits, respectively. Observe that the moments depend on the starting point as well 
as on the initial direction of the neutron, and level off to an asymptotic value, which is 
the signature of the process going to ultimate extinction. Furthermore, again for (£y )t 
we have {iv)t~ — {^v)^ to min(r+,r"). 

Observe finally that the rod model has been widely adopted to describe, among 
others, gas dynamics or biological species migration [I5l[l6],[171[l8l[l9], so that the results 
described in this Section concerning an application in reactor physics could perhaps be 
of interest in other areas of science as well. 
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7. Summary and conclusions 

In this paper we liave examined tlie statistics of travelled lengths and number of collisions 
for branching exponential flights. Moment formulas have been derived by resorting 
to the backward Feynman-Kac formalism, based on a minimal number of simplifying 
hypotheses. Moreover, we have shown that this same formalism can be extended with 
slight modifications to the analysis of other physical observables, such as the survival 
probability. The proposed formulas have been compared to Monte Carlo simulations for 
an example of one-dimensional transport inspired by reactor physics, and an excellent 
agreement was found. A generalization to more complex transport problems would also 
be possible, by relaxing for instance the requirement on the isotropy of the scattering 
kernel [50] and by introducing energy and space dependent cross sections. We conclude 
by observing that, although throughout this paper we have used our knowledge of ip 
and (p so to assess the statistical properties of the observables iy and ny, the inverse is 
also possible: for instance, knowledge of the average total travelled length in V allows 
inferring the average particle flux in that volume, and similarly knowledge of the total 
number of visits to V allow inferring the average collision density in that volume. 

Appendix A. The Feynman-Kac equations for Qt and Wt 

Consider a single walker initially at fq, ujq at observation time t = 0. We want to write 
an equation for the moment generating function Qt = {e~^^^^*^). Assume an observation 
time t + dt: this can be split into a first interval, from to dt, and then a second interval 
from dt to t + dt. The only requirement is that the process is Markovian: after dt the 
particle continues its path without memory of the past. We start by observing that in 
a vanishing small time interval dt, from the definition of the underlying process, only 
two mutually exclusive events are possible: either the particle does not interact with 
the medium, in which case the walker keeps going in the same direction by a space 
interval dr^ = vuJodt, or the particle interacts, in which case the walker disappears and 
gives rise to k descendants at the same position, with random directions uj'^ obeying 
the same probability density C. From the definition of the cross section a, the former 
event happens at a rate 1 — avdt, whereas the latter at a rate avdt. If no particles are 
emitted (po), the trajectory is terminated, and no further contribution is added to iy 
When k > 1 (identical) particles are generated, the probability that the contribution to 
the total travelled length coming from each walker adds up precisely to iy is given by 
the convolution of the probability that the first particle spends a length £i, the second 
£2, and the k-th a length iy — £1 — £2 — ■ ■ ■■ In the transformed space, the convolution 
products amount to a simple product of generating functions. This simple argument 
leads to the equation 

Qt+dt{s\ro,oJo) = (1 - (Ti;rft)e-^^^('-°''^°)'^*gi(s|ro + rfro,a;o) + 

avdt [po+Pi{Qtis\ro,u;[)) + p2(Qt(s|ro, a;2i)Qt(s|ro, a;22)) H ] , (A.l) 
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where brackets denote expectation with respect to the random directions cj'^. Now, if we 
suppose that the descendant directions are independent, the expectation of a product 
of random variables becomes the product of expectations, so that we are led to 

where G[z\ = Pq + piz + p2Z^ + ■ ■ ■ is the generating function associated to pk- Observe 
that the average over the random directions can be expressed in terms of the associated 
probability density as 

((5t(s|ro,a;o)) = j C*{u;'q u;o\ro)Qtis\ro, u;'Q)duj'Q, (A.3) 

where C* is formally the adjoint density with respect to C. As a shorthand, we will 
denote C*{Qt} = J C*{uj'q — )■ uJo\TQ)Qt{s\rQ,uj'Q)duj'Q. Now, when dt is small, at the 
leading order we have Qt{s\ro + dTQ,u)o) = Qt{s\rQ,ujQ) + vujq ■ VroQtdt + ■ ■ ■, along 
the direction of loq. Furthermore, for vanishing dt we have exp{—svV{ro,ujQ)dt) = 
1 — svV{rQ, ujQ)dt + ■ ■ ■. By recollecting all terms we then get 

Qt+dt = Qt+vu)Q-Wr,Qtdt-vaQtdt-svV{vQ,ijjQ)dt+avdtG [C*{Qt}] .(A.4) 

By dividing by vdt and taking the limit rft — )■ 0, we finally obtain the backward Feynman- 
Kac equation for the moment generating function Qt, namely 

^-^Qt = ivo ■ Vr,Qt - oQ^ - sV(ro,a;o)gt + aG[C*{gi}]. (A.5) 

Now, from the same argument as above, it follows that the probability generating 
function Wt = (s™^(*)) satisfies 

Wt+dt{s\ro, ujq) = (1 - (rvdt)Wt{s\ro + dro, ujq) + 

avdt[po + Pi{Wt{s\ro,uj[)) + p2{Wt{s\ro,u:'^^)Wt{s\ro,u:'^2)) + ■■■], (A.6) 

hence 

^^Wt = cjo ■ VroWt - aWt + aG[C*{Wt}]. (A.7) 
Appendix B. The discrete Feynman-Kac equation for Qn 

Let a single walker be initially at tq, (jJq. When considering generations, it is more 
convenient to begin with a single particle entering its first collision with coordinates 
ri,a;o, and we denote by Qn{u\Yi^ujQ) = (e^""^*^"'^)(ri, cjq) the corresponding moment 
generating function. The separation between ri and tq, at a first glance somewhat 
artificial, is actually due to the special role of the source: a particle emitted from the 
source is just transported to the first collision point, and can not be absorbed nor 
multiplied at tq [lOlHT]. When entering the collision at ri, k particles are created, with 
probability pk, and random directions. Exponential flights are Markovian at collision 
points, which allows splitting each of the k subsequent trajectories into a first jump, 
from ri to r'^ in direction uj'^^, and a branching path from r'^ to the positions held at the 
(n + l)-th generation. The displacement r'^ — ri obeys the jump length density T, and 
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the direction uj'^ the density C. li k = 0, the trajectory ends at ri and there will be no 
further events contributing to ny. Hence, we have 

g„,+i(^i|ri,a;o) =Poe-"^('-^''^«) +Pie-"^('-^''^°)(gn(n|r;,a;;)) + 

where expectation is taken with respect to the random displacements and directions, 
and the term e"'"^*^'"^''^") can be singled out because it is not stochastic. The terms at 
the right hand side in Eq. IB.ll can be understood as follows: the probability that k 
identical and indistinguishable particles (born at ri) give rise to ny collisions in V is 
given by the convolution product that the first makes ni collisions, the second n2, 
and the fc-th ny — ni — n2 — ■ ■ ■■ In the transformed space, this convolution becomes 
a simple product of generating functions. If we assume that the descendant particles 
are independent, the expectation of the products in Eq. IB. II becomes the product of the 
expectations. Now, observe that the average over displacements and directions can be 
expressed in terms of the associated densities, namely 

{Qniu\r'^,u;'^)) = j C*{u}'q ^ u;o\ri) J T*{v[ ^ ri\u;'f))Qniu\v[,uj'Q)dv[du;'Q, iB.2) 

where T* is the adjoint density associated to T j40j. Intuitively, T* displaces the walker 
backward in time. Observe in particular that the first collision coordinates ti,ujq obey 
the probability density T*(ri — )■ rol^o), namely, 

Qn{u\ro,u;o) = j Qniu\ri,uJo)T*{ri ro\ujo)dri. (B.3) 

Therefore, by using C*{Qn} = f C*i<^o — ^ ^o\''''i)Qn{u\ri,^'o)duj'Q as above, we obtain 
the discrete Feynman-Kac equation in integral form, namely 

Qn+iiu\v,, ivo) = e-«^(^^''^»)G' [C*{Q„}] . (B.4) 

Finally, by integrating over T* both sides of Eq. IB. 41 we get 

g„,+i(^i|ro,a;o) = j dr^T*{r^ ^ roMe-''''^'''''^'^G[C*{Qn}] . (B.5) 

It can be shown [23] that integral equations in the form 



/(ro,Wo) = J driT*{ri Fq 1 0^0)5' (ri, a;o) (B.6) 
can be equivalently recast into 

'^^o ■ Vro/(ro, ujo) - o-/(ro, cjq) + (yg{ro, ujq) = 0. (B.7) 
Therefore, Eq. IB. 41 gives the discrete Feynman-Kac equation in integro-differential form 

-^0- VroQn+i(^|ro, u:o)+aQn+i{u\ro, u:o) = ae-^'^'^'^'^'^^G [C*{Q„}] .(B.8) 

Now, from the same argument as above, it follows that the probability generating 
function Wn = (■u™^^"')) satisfies 

))+---,(B 

with iy„(ri,a;o) = (M'"^("))(ri,Lc;o). Finally, 

- Wo ■ Vrol^n+i(M|ro, a;o) + aiy„+i(M|ro, = aC [C*{Wn}] . (B.IO) 
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Appendix C. The rod model equations 



From Eq. [161 the average traveUed length obeys 
ld_ 

V dt 

which gives then 



-{e],)t = c*{ey)t + v{xo){i'v)t, 



(C.l) 



ld_ 

V dt 



= -^/v)t - + ^ mt + + ^(^o). (C.2) 



For the second moment of the travelled length we have 
1 d 



{fy)t = C*{fy)t + 2V{xo){fv}t + az/2C*{(4)t}^ 



which gives 



V dt 

V at 



(C.3) 



a 



d 



+2V{xo){fy)t + Y i^^v^t + {^v)iy 

+2ViXo){fy), 



dxo 



From Eq. 
which gives 



* ' 4 

the average number of visits obeys 

- '^^O • Vro(ny)„+l + a{ny)n+l = CTl^lC* {{Uy) n} + V{Xo), 



(C.4) 
(C.5) 



d 
dxo 



{K)t+{<)n)+V{Xo), 



For the second moment of the number of visits we have 

+2az/iy(xo)C*{(nt.)4 + aV{x^), 
which leads to 



(C.6) 



(C.7) 



d 



+ ^ {{n'v)t + (On)' + ^ + {<)-) + crVixo) 

+^ {{<)i + iOnY + ^ («);[ + «);) + ^^(^o). 



(Ci 
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